This report presents an overview of local and statewide Housing Choice Voucher (HCV) usage to better understand where source of income discrimination may be occurring.
This report was built using open-source programming, tools and datasets.
All datasets used in this report are publicly-available and from administrative sources.
# Import Housing Choice Voucher Data
hcv <- read.csv("data/HUDPicture_2017_HCV.csv", stringsAsFactors = FALSE) %>%
filter(Id2 != "42XXX999999") %>%
mutate(tract = Id2) %>%
mutate(county = str_sub(tract, 1, 5)) %>%
select(tract, county, hcv_sub_units)
tenure <- read.csv("data/ACS_17_5YR_B25003_with_ann_TENURE.csv", stringsAsFactors = FALSE, colClasses=c("Id2"="character")) %>%
mutate(tract = Id2) %>%
select(-Id, -Id2)
hcv <- hcv %>%
left_join(tenure, by = "tract") %>%
mutate(hcv_rate = ifelse(tothh < 10, NA, 100 * hcv_sub_units/tothh))
hcv.county <- hcv %>%
group_by(county) %>%
summarise(hcv_sub_units = sum(hcv_sub_units, na.rm = TRUE),
tothh = sum(tothh, na.rm = TRUE)) %>%
as.data.frame() %>%
mutate(hcv_rate = 100 * hcv_sub_units/tothh)
# Import Census Data
hisp <- read.csv("data/ACS_17_5YR_B03002_with_ann_HISP.csv", stringsAsFactors = FALSE)
pov <- read.csv("data/ACS_17_5YR_S1701_with_ann_POV.csv", stringsAsFactors = FALSE)
hisp <- hisp %>%
mutate(tract = as.character(Id2)) %>%
mutate(tract_poc = totpop - nhisp_wh) %>%
mutate(tract_pct_poc = 100*(totpop - nhisp_wh)/totpop) %>%
select(tract, totpop, tract_poc, tract_pct_poc)
pov <- pov %>%
mutate(tract = as.character(Id2)) %>%
mutate(tract_pct_pov = 100 * (belpov/totpovstatus)) %>%
select(tract, totpovstatus, belpov, tract_pct_pov)
# Import Tract Shapefile for Philadelphia
tracts <- st_read("shps/tl_2018_42_tract.shp", layer = "tl_2018_42_tract", stringsAsFactors = FALSE) %>%
mutate(tract = GEOID) %>%
select(tract, geometry)
# Import County Shapefile
counties <- st_read("shps/PA_Counties.shp", layer = "PA_Counties", stringsAsFactors = FALSE) %>%
mutate(state = str_sub(GEOID, 1, 2)) %>%
filter(state == "42") %>%
mutate(county = GEOID) %>%
select(county, geometry, NAME)
# Join data to tracts and counties
tracts <- left_join(tracts, hcv, by = "tract")
tracts[, 3:7][is.na(tracts[, 3:7])] <- 0 # Replace NAs with 0s
tracts <- tracts %>%
mutate(county = str_sub(tract, 1, 5)) %>%
left_join(hisp, by = "tract") %>%
left_join(pov, by = "tract") %>%
mutate(recap = ifelse(tract_pct_poc>=50 & tract_pct_pov >= 40, "RECAP", "NOT RECAP")) # RECAP definition by HUD: https://data.world/hud/recap
counties <- left_join(counties, hcv.county, by = "county")
# Import HUD Small Area FMRs and Zip Codes
safmr <- read.csv("data/hud_phila_small_fmrs_fy2019.csv", stringsAsFactors = FALSE)
fmr <- read.csv("data/hud_phila_fmrs_fy2015.csv", stringsAsFactors = FALSE) #https://www.huduser.gov/portal/datasets/fmr/fmrs/FY2015_code/2015summary.odn
zcta <- st_read("shps/Phila_ZCTA_WGS84.shp", layer = "Phila_ZCTA_WGS84", stringsAsFactors = FALSE)
# Import Median Gross Rent Data
rent <- read.csv("data/ACS_17_5YR_B25064_with_ann_RENT.csv", stringsAsFactors = FALSE) %>%
mutate(high_moe = ifelse(moe_median_rent/median_rent > 0.40, 1, 0)) %>%
mutate(tract = as.character(GEOID))This section uses tract-level data from HUD on voucher holders in 2017 to show the spatial distribution.
The following is a dot density map of voucher holders across Pennsylvania. Each dot represents 10 households with vouchers and is randomly placed within it’s respective Census tract (the dots do not represent the exact location of households). Click on each county to see the number and percentage of voucher holders.
dots <- suppressMessages(st_sample(tracts, size = round(tracts$hcv_sub_units/10), type = "random")) # each dot = 10 hcv units
popup <- paste0("County: ", counties$NAME, "<br>", "Voucher Rate: ", round(counties$hcv_rate, 2), "<br>", "Vouchers: ", counties$hcv_sub_units)
m <- leaflet(options = leafletOptions(preferCanvas = TRUE)) %>%
addProviderTiles("Esri.WorldGrayCanvas", options = providerTileOptions(
updateWhenZooming = FALSE, # map won't update tiles until zoom is done
updateWhenIdle = TRUE # map won't load new tiles when panning
)) %>%
addCircleMarkers(data = dots, weight = 0, color = "#18BC9B", radius = 3) %>%
addPolygons(data = counties, fill = FALSE, stroke = TRUE, color = "#626468", weight = 1, popup = ~popup) %>%
addLegend("bottomright", colors= "#18BC9B", labels="1 dot = 10 HCVs", title="Vouchers in PA Counties")
mThe following map shows the estimated percent of voucher holders by neighborhood in Philadelphia. Click on individual neighborhoods to see the number and percent of voucher holders in that area. Neighborhood estimates were produced from tract-level estimates by applying area-weights.
# Calculate HCVs by Neighborhood
tract.neigh <- read.dbf("shps/Tracts_Neigh_Int.dbf", as.is = TRUE) %>%
mutate(tract = GEOID) %>%
mutate(neigh = MAPNAME) %>%
mutate(int_fract = AreaInt_ft/Area_ft) %>%
select(tract, neigh, Area_ft, AreaInt_ft, int_fract) %>%
left_join(hcv, by = "tract") %>%
mutate(hcv_weighted = hcv_sub_units * int_fract) %>%
mutate(hh_weighted = tothh * int_fract)
neigh.hcv <- tract.neigh %>%
group_by(neigh) %>%
summarise(
hcv = sum(hcv_weighted, na.rm = TRUE),
hh = sum(hh_weighted, na.rm = TRUE)
) %>%
as.data.frame() %>%
mutate(hcv_rate = ifelse(hh < 10, NA, 100*hcv/hh))
neigh <- st_read("shps/Neighborhoods_Phila_WGS84.shp", layer = "Neighborhoods_Phila_WGS84", stringsAsFactors = FALSE) %>%
mutate(neigh = MAPNAME) %>%
select(neigh, geometry) %>%
left_join(neigh.hcv, by = "neigh")## Reading layer `Neighborhoods_Phila_WGS84' from data source `T:\GitHub\soi\shps\Neighborhoods_Phila_WGS84.shp' using driver `ESRI Shapefile'
## Simple feature collection with 158 features and 5 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: -75.28027 ymin: 39.86701 xmax: -74.95576 ymax: 40.138
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defsbins <- quantile(neigh$hcv_rate, seq(0,1, by = 1/5), na.rm = TRUE)
pal <- colorBin("YlOrRd", domain = neigh$hcv_rate, bins = bins)
popup <- paste0("Neighborhood: ", neigh$neigh, "<br>", "Voucher Rate: ", round(neigh$hcv_rate, 2), "<br>", "Vouchers: ", round(neigh$hcv, 0))
m2 <- leaflet(neigh) %>%
setView(lng = -75.140012, lat = 39.999786, zoom = 11) %>%
addProviderTiles("Esri.WorldGrayCanvas") %>% #Stamen.TonerBackground
addPolygons(color = "#ffffff", weight = 1, smoothFactor = 0.5, opacity = 0.8, fillOpacity = 0.5,
fillColor = ~pal(hcv_rate), popup = ~popup,
highlightOptions = highlightOptions(color = "#444444", weight = 2, bringToFront = TRUE)) %>%
addLegend(pal = pal, values = ~hcv_rate, opacity = 0.7, title = "Voucher Holders </br> per 100 Households", position = "bottomright")
m2neigh.tbl <- neigh
st_geometry(neigh.tbl) <- NULL
neigh.tbl <- neigh.tbl %>%
mutate(hcv = round(hcv),
hh = round(hh),
hcv_rate = round(hcv_rate,2))
datatable(neigh.tbl, rownames = FALSE, colnames = c("Neighborhood", "HCV", "Households", "HCV per 100 Households"))HUD established the definition of Racially and Ethnically Concentrated Areas of Poverty (RECAP) as Census tracts with:
See: https://data.world/hud/recap for more information.
This table presents the number and percent of tenants with vouchers in RECAP tracts vs. non-RECAP tracts.
recap <- tracts %>%
group_by(county, recap) %>%
summarise(
total_households = sum(tothh, na.rm = TRUE),
total_hcv = sum(hcv_sub_units, na.rm = TRUE),
ave_hcv_rate = mean(hcv_rate, na.rm = TRUE),
ave_pct_poc = mean(tract_pct_poc, na.rm = TRUE)
) %>%
as.data.frame() %>%
left_join(counties, by = "county") %>%
filter(is.na(recap)==FALSE) %>%
select(county, NAME, recap, total_hcv) %>%
spread(recap, total_hcv) %>% # convert to wide
as.data.frame() %>%
rename(not_recap = `NOT RECAP`) %>%
rename(recap = RECAP) %>%
mutate(pct_recap = round(100*(recap/(recap+not_recap)),2))
datatable(recap, rownames = FALSE, colnames = c("County ID", "County", "HCVs not in RECAPs", "HCVs in RECAPs", "% HCVs in RECAPs"))Results for Philadelphia:
phl.tracts <- tracts %>%
filter(county == "42101")
phl.recap <- phl.tracts %>%
filter(recap == "RECAP")
bins <- quantile(phl.tracts$hcv_rate, seq(0,1, by = 1/5), na.rm = TRUE)
pal <- colorBin("YlOrRd", domain = phl.tracts$hcv_rate, bins = bins)
popup <- paste0("Voucher Rate: ", round(phl.tracts$hcv_rate, 2), "<br>", "Vouchers: ", round(phl.tracts$hcv_sub_units, 0))
m3 <- leaflet(phl.tracts) %>%
setView(lng = -75.140012, lat = 39.999786, zoom = 11) %>%
addProviderTiles("Esri.WorldGrayCanvas", group = "Base") %>%
addPolygons(group = "Voucher Rate",
color = "#ffffff",
weight = 1, opacity = 0.5,
fillOpacity = 0.5, fillColor = ~pal(hcv_rate),
popup = ~popup,
highlightOptions = highlightOptions(color = "#444444", weight = 2)) %>%
addPolygons(data = phl.recap,
group = "RECAP",
color = "#3b485e", weight = 3,
fill = "#ffffff", fillOpacity = 0.2) %>%
addLegend(pal = pal, values = ~hcv_rate, opacity = 0.7, title = "Voucher Holders </br> per 100 Households", position = "bottomright") %>%
addLegend(data = phl.recap, colors = "#3b485e", labels = "RECAP Tracts", position = "bottomright") %>%
addLayersControl(
baseGroups = "Base",
overlayGroups = c("Voucher Rate", "RECAP"),
options = layersControlOptions(collapsed = FALSE)
)
m3Based on NYC analysis by the Commision on Human Rights and Mayor’s Office of Data Analytics.
NYC Methodology:
In late 2016, HUD began using small area Fair Market Rents (FMRs) in certain metropolitan areas to set payment standards for Housing Choice Voucher recipients. In the Philadelphia-Camden-Wilmington metro area, zip codes are used as small area FMRs.
Small area FMRs were an attempt to correct for the error in using county or metro-wide Fair Market Rent estimates. In theory, small area FMRs would be closer to the on-the-ground rent reality; however, zip codes are still large geographies that do not match perfectly to more granular neighborhood-level rental prices. For example, in gentrifying areas, small area FMRs may underestimate rents, further disincentivizing landlords from renting to tenants with vouchers.
The following map shows the difference between 2019 small area FMRs for a 2-bedroom unit and the 2015 metro-wide FMR for a 2-bedroom unit ($1,156). A positive number reflects areas
safmr <- safmr %>%
mutate(efficiency = gsub("\\$", "", efficiency)) %>%
mutate(br_1 = gsub("\\$", "", br_1)) %>%
mutate(br_2 = gsub("\\$", "", br_2)) %>%
mutate(br_3 = gsub("\\$", "", br_3)) %>%
mutate(br_4 = gsub("\\$", "", br_4)) %>%
mutate(efficiency = gsub(",", "", efficiency)) %>%
mutate(br_1 = gsub(",", "", br_1)) %>%
mutate(br_2 = gsub(",", "", br_2)) %>%
mutate(br_3 = gsub(",", "", br_3)) %>%
mutate(br_4 = gsub(",", "", br_4)) %>%
mutate(efficiency = as.integer(efficiency)) %>%
mutate(br_1 = as.integer(br_1)) %>%
mutate(br_2 = as.integer(br_2)) %>%
mutate(br_3 = as.integer(br_3)) %>%
mutate(br_4 = as.integer(br_4)) %>%
mutate(br2_diff = br_2 - 1156)
zcta.fmr <- zcta %>%
mutate(zip = as.integer(GEOID10)) %>%
select(zip, geometry) %>%
left_join(safmr, by = "zip")
bins <- c(-166, -75, -5, 5, 200, 644)
pal <- colorBin("PuOr", domain = zcta.fmr$br2_diff, bins = bins)
popup <- paste0("Difference $: ", zcta.fmr$br2_diff)
m4 <- leaflet(zcta.fmr) %>%
setView(lng = -75.140012, lat = 39.999786, zoom = 11) %>%
addProviderTiles("Esri.WorldGrayCanvas") %>%
addPolygons(color = "#ffffff", weight = 1, smoothFactor = 0.5, opacity = 0.8, fillOpacity = 0.5,
fillColor = ~pal(br2_diff), popup = ~popup,
highlightOptions = highlightOptions(color = "#444444", weight = 2)) %>%
addLegend(pal = pal, values = ~br2_diff, opacity = 0.7, title = "FMR Difference", position = "bottomright")## Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
## Need '+proj=longlat +datum=WGS84'm4# lit <- read.csv("SOI Literature Summary.csv", stringsAsFactors = FALSE)
# datatable(lit, rownames = FALSE)Methodology:
phl.tracts.rent <- phl.tracts %>%
left_join(rent, by = "tract") %>%
filter(high_moe != 1)
bins <- quantile(phl.tracts.rent$median_rent, seq(0,1, by = 1/5), na.rm = TRUE)
pal <- colorBin("YlOrRd", domain = phl.tracts.rent$median_rent, bins = bins)
popup <- paste0("Rent: ", phl.tracts.rent$median_rent)
m6 <- leaflet(phl.tracts.rent) %>%
setView(lng = -75.140012, lat = 39.999786, zoom = 11) %>%
addProviderTiles("Esri.WorldGrayCanvas") %>%
addPolygons(color = "#ffffff", weight = 1, smoothFactor = 0.5, opacity = 0.8, fillOpacity = 0.5,
fillColor = ~pal(median_rent), popup = ~popup,
highlightOptions = highlightOptions(color = "#444444", weight = 2)) %>%
addLegend(pal = pal, values = ~median_rent, opacity = 0.7, title = "Median Rent by Tract", position = "bottomright")## Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
## Need '+proj=longlat +datum=WGS84'm6